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We consider the formation of defects in a non equilibrium second order phase transition induced 
by an instantaneous quench to zero temperature in a type II superconductor. We perform a full 
non-linear simulation where we follow the evolution in time of the local order parameter field. We 
determine how far into the phase transition theoretical estimates of the defect density based on 
the gaussian approximation yield a reliable prediction for the actual density. We also characterize 
quantitatively some aspects of the out of equilibrium phase transition. 
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I. INTRODUCTION 



X> 

, The objective of this paper is to study the formation of defects in a non equilibrium second order phase transition by 
means of a numerical solution of the full dynamical equations, and to compare the results with theoretical predictions 
I/") to be found in the literature 

O^I ■ Topological defects are a common occurrence in symmetry broken field theories [0 , with far reaching consequences 
in condensed matter physics || , particle physics ||[7J and cosmology || . The equilibrium structure of defects is deeply 
rooted in the topological aspects of the theory and it is well understood. The dynamical formation of defects in the 
process of a non equilibrium phase transition, on the other hand, only recently has been the subject of a systematic 
\ analysis. 

£f) • The issue at stake is how many defects are to be formed, as a function of both the dynamics of the system and the 
macroscopic parameters characterizing the transition, such as cooling rates. A simple, order of magnitude estimate 
■ is based on the observation that once the system is cold enough, defects will be unable to either form or disappear 
through thermal activation, so they will be essentially frozen into existence. This happens at the so-called Ginzburg 
temperature, and leads to the prediction that the typical distance between defects is of the order of magnitude of the 
' correlation length at this temperature || . 
Oh, This picture of defect dynamics has been criticized by Zurek as downplaying the non equilibrium features of the 
process Jl(],|ll[]. According to this author, the freezing of defects occurs at a much higher temperature, when the 
relaxation time of the system becomes large and the system effectively decouples from its environment. Zurek's 
1 arguments usually lead to higher defect densities than previously though, a prediction confirmed by some experiments 
[ p2| , p^[ (see however |l4[). They are supported also by numerical simulations in one and two dimensional Ginzburg- 
Landau systems with broken global symmetries . 

Complementary to the search for a qualitative understanding of defect formation, some authors have attempted to 
derive the density of defects from a first principles description of the dynamics. As a matter of fact, there is a rigorous 
result linking the density of defects to the equal time correlation function for the order parameter, valid whenever 
the field probability distribution is Gaussian QJ^]. Under the Gaussian approximation, then, finding the density of 
defects is reduced to solving the dynamics for the correlation function ||. This approach, which has been recently 
used to show the validity of Zurek's estimates in a quantum field theoretic model , is usually though to be correct 
early in the development of the phase transition p7| . 

In order to evaluate how reliable this kind of argument really is, therefore, we must know how far in the phase 
transition the Gaussian approximation may be trusted. Since the development of the phase transition is an essentially 
nonlinear process, we must expect that non Gaussian correlations will be created by the dynamics itself, even if 
suppressed in the initial conditions. There will be a competition between the characteristic growth time due to the 
spinodal instability Jl^,|l^l , and a variety of dynamical times describing the building of correlations through non linear 
interaction of fluctuations; however, the usual Gaussian models do not give us a clue about what the latter might be. 
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The objective of this paper is to give a tentative answer to the question of the reliability of defect density estimates 
based on the Gaussian approximation, by presenting a fully nonlinear simulation of a phase transition where we have 
measured the Gaussianity of the order parameter, the Gaussian prediction for the density of defects, and the actual 
density, as functions of time. Of course, as long as the held is actually Gaussian, the rigorous analytic prediction for 
the defect density is validated. The nontrivial issue is how long the held probability density remains Gaussian, and 
whether the Gaussian prediction continues to hold beyond that point. Concretely, we find that even when gaussianity 
of the ensemble ceases to hold, the prediction remains valid, to break down about the time defects are definitely 
formed. 

Our simulations follow the unfolding of an instantaneous quench to zero temperature in a two dimensional type 
II superconductor, described by the time-dependent Landau-Ginzburg equations derived by Gor'kov and Eliashberg 
pp| (sec also ]2l|-p4|| ) . This model is a scalar field theory interacting with a U(l) gauge field (see |]ll|) but with the 
presence of a normal current, besides the supercurrent, and a first order (rather than second one) equation for the 
vector potential A (see below). 

In order to perform the simulations, we have discrctized the model and placed it on a square lattice, with periodic 
boundary conditions. Our discrete model still has gauge symmetry, which is preserved by the evolution. We have 
further considered a variety of lattice sizes and initial conditions, thus making sure that our results truly reflect the 
physics of the system. We have identified defects by measuring the winding number of the field around each lattice 
plaquette, thus avoiding the uncertainties related to identifying defects from zeroes in the order parameter p5(| . 

We find that the evolution of a typical quench goes through three well defined regimes. The early development is 
dominated by the exponential growth of the order parameter; the different modes of the field evolve independently, 
and the field remains Gaussian. In this regime, the Halperin - Mazenko - Liu (HML) prediction for the defect density 
is very accurate. This regime ends when the order parameter reaches about a tenth of its equilibrium value. 

The second regime is a transitional epoch dominated by the actual formation of the defects. During this transitional 
epoch, the field departs from gaussianity in a significative way, but the HML prediction is still a good approximation 
to the actual density. Finally, in the late time regime both the Gaussian approximation for the order parameter and 
the HML prediction are unreliable. 

We may therefore conclude that, as argued by Karra and Rivers j|, the HML prediction holds early in the devel- 
opment of the phase transition, being a very accurate estimate of the actual density until the time the defects may be 
considered as definitely formed, which is also the time when the order parameter reaches about a half of its equilibrium 
value. It is therefore a suitable means to estimate the initial conditions for the subsequent evolution of the defect 
network, as determined by defect - defect interactions and changes in the environment, such as the expansion of the 
Universe in cosmological applications j^] . Our results validate previous analysis of non equilibrium defect formation, 



The paper is organized as follows. The first part introduces the theoretical prediction of the defect density for 
a gaussian quench. The second part shows the time-dependent Ginzburg-Landau model, and the numerical details 
related to its implementation. The third part describes the resulting quenches and the conclusions, and the paper 
ends with some final remarks. 



At a qualitative level, the formation of topological defects is well understood through the symmetry breaking 
mechanism. For a complex scalar field, the true ground-state manifold of the field is SI and the phase of the field at 
different points can be different. If eventually the winding number of the field along a closed loop is not zero, then 
topological defects are trapped inside. Because of topological considerations, isolated defects at low temperature are 
stable, although defects may interact with each other and annihilate, migrate to the boundaries of the system, or, in 
cosmology, decay through gravitational radiation. 

We are searching for the monopoles of one complex order parameter field $ in two dimensions. This will be identified 
with the zeroes of the field with non-trivial winding number. If they are located at xi, X2, X3, we obtain for the 
total and topological densities 



such as I |JT§|. 



II. THEORETICAL PREDICTION 




(1) 




Where n.; is the winding number of each defect, i.e. its topological charge. 
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The total density of defects is obtained through the relation between zeros of the field and the field itself, i.e. the 
Jacobian ||, 



i(t) = (p(x)) = J V<Pp t [$] 5 2 [$] |e jfe a j $i(x)9 fc $ 2 (x)| 



(2) 



with e 12 = —£2i = 1 (otherwise zero) and $ = (<&i + «$2)/v / 2, Pt being the probability density of the different field 
configurations. 

For the Gaussian model we have ($ a (x)) = = ($ a (x)<9,-$&(x)). Assuming also that the equal-time Wightman 
function (<j> (x)<I>;,(y)) = VF a h(|x — y| ;t) = 5 a t,W (|x — y| ; t) is the only non-vanishing correlation function and it is 
diagonal, then 



n(t) = r-(-/"(0;t)) 

Z7T 



where 



f(r;t) = 



W(r;t) 
W(0;t) 



(3) 



(4) 



and derivatives are taken with respect to r. 

To compute the r.h.s. of (|J) we consider the Fourier transform of the field, $(k) = J <ixexp(ik.x)<I>(x) in order to 
obtain 



f"(0;t) 



fdkk 2 


$(k) 2 


fdk 


i(k) 2 



(5) 



Our goal is to test the range of validity of the relation (0) by measuring both sides of this equation independently. 



III. THE MODEL 



A. Theory 



The time-dependent Ginzburg-Landau equation describes the time-space dependence of the order parameter of a 
superconductor pOf . Normalized in the form adopted by Hu and Thompson, it reads {2lLp3| 
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J = o- 
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A-/(r,t)=0 



he 2 



47rA^ F 



(6) 

(7) 
(8) 



With the Maxwell equations coupling the electromagnetic potentials to charge and current densities, they provide 
the full set of evolution equations. Here, D is the normal-state diffusion constant, a is the normal-state conductivity 

given by a — 48 ^2 p and j are the charge and current densities. / is a finite temperature random driving force; 
since we shall consider a quench to zero temperature, it will eventually be set to zero (see below). Furthermore, the 
order parameter is divided by its equilibrium value Aqo = nk^2 (T c 2 — T 2 ), where T c is the critical temperature. Note 
that this is a temperature-dependent parametrization. 

A and if are the vector and scalar potentials respectively, and tp is the electrochemical potential divided by the 
electronic charge. Assuming p — 0, that is, absence of net charge at the grid scale, results in t/j = cp (see ref. 

This set of equations is invariant under the gauge transformation 
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A 



A V X 
1 9x 

c flt 



(9) 



Aexp 



.2e 



From the microscopic theory we have the relationships 



47rA(T) 2 (7 £(T) 2 



tt/j 



where £(T) = £(0) 



1 



-1/2 



12D 96/c B T c 



T 



L GLQ 
12 



T 



(10) 



and A(T) = A(0) 



1 



-1/2 



length and magnetic penetration depth respectively, and iGLO = 8fc T 
uniform mode at zero temperature. 

We can write our model in terms of dimensionless variables as follows j24 

ttS, taw 



are the temperature dependent correlation (coherence) 
is the characteristic relaxation time of the 



t — ► t to with to 
r -> r £(0) 



96fc s T c 12 
he 



(11) 



^4 ^ ; 7^ with $ 

2tt£ (0) 2e 



27TCt 

c$ 



to obtain 



8^ 2 £(0) 

/-/£(o) 2 



^A + = -i [(iv + A) 2 A + (1 - T) (|A| 2 - l) A - / 
— A + Vtp = (1 - T) Re [A* (-iV - A) A] - k 2 V x (V x A) 



(12) 



where k = A(T)/£(T) is the temperature independent Ginzburg-Landau parameter which characterizes the supercon- 
ductor. For a type II superconductor we have k > l/\/2. We have chosen k = \[2. 
The gauge freedom allows us to choose tp = 0. 

Since we are interested in instantaneous quenches towards zero temperature, we set T = 0, and / = f p6[ . 
Furthermore we must prepare the system in some thermal equilibrium configuration. This means to generate a set of 
initial conditions corresponding to a thermal distribution of modes u% '■ 



A(k) 



(A(k)A(O) 





1 2m* 



(13) 



V h 2 fc 2 + l/£(T) 5 



with a cutoff when k sa £ 1 beyond which GL theory is not valid. Here m* = 2m e , is the mass of the coupled electrons. 
In terms of dimensionless variables, 



A(k) 



fi T 



with /i = 



2k B T c m*i (0) 



V fc 2 + £~ 2 ' fr« 
The factor \i is clearly substance dependent and will be chosen later on. 



(14) 
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B. Implementation 



The discrete version of the gauge transformation is obtained through 



Al _» A3 _ 



(15) 



A j -> A 3 exp [-ix°] 



where p, stands for a direction and j for a site in the lattice. In order to obtain a discrete version of (|l2]), invariant 
under ( p"5| ) we employ the usual link variables technique from Lattice QCD |?j as in J24|] : 



f/; ir2 = exp 



1-2 

-i J A^dp 

ri 



discretized U 3 ' ■ >+l1 — exp [— iA^a^] 



(16) 



The differential operators become 



1 d 



i dx^ 

i d 

i dx u 



A, 



A, 



Ijj+H, 3/±3+H _ /±j 



3^+^ - 2A j + U?r»' A J ' _/i 



and the finite difference equations to solve are ( A° — ffle %e ): 

Re A J =— (p J+x cosf- A%a x + - 2 Re A 
12 [ v 

-(l-T)(p j2 -l)ReA J ] 
for the real part of the order parameter field, and 



' J- - m>;,l< •'•„, + W ')) — 



(17) 
(18) 

(19) 



:C = (1 - T)pV+* sm(-a x Ai - W + 6 j+x ) — 



(20) 



A{+*+y - - A{ +x + A( Al+y - 2Ai + At* 



for the x component of the gauge field. We choose here to work directly with the fields, but alternatively the link 
variables can be used j p4j . 

We evolve this equation with a simple Euler scheme, taking time-steps empirically chosen to be h = io/128 (t 
being the time scale defined in eq. |TT| ) , and imposing periodic boundary conditions. The choice of the time step is 
very critical because of the very different and variable time scales involved in this kind of simulation. 

We tested grids of N 2 sites, with N = 128. 25G and 512. It is convenient to employ larger grids, not only because of 
less granularity in the observed density, but also because it is possible to achieve sufficient statistics with fewer runs 
for ensemble. We have made about 20 runs for each ensemble, which means that the dispersions of field and defect 
density are in the ~ 3% level. 

We choose the net parameters a x — a y — £o/2, in order to resolve adequately the shape of the defects, which are 
expected to have a final size d ~ £o- 

The initial conditions were set in two different ways. 

We obtain a thermal distribution of modes, employing expression (|l4|) as the dispersion of a gaussian distribution 
of the mode amplitudes, with £ = £o/2, that is a temperature equal to five times T c . The cutoff was set up at the 
maximum radius in fc-space, i.e. k max = The results are cutoff independent, in any case, due to the rapid decay 
of the short wavelength modes in the first steps of the quench. At the end of the quench, these modes grow again in 
order to define the final shape of the defects. Following |24], we tested p = 10~ 2 , 10~ 3 , 10~ 4 and 10~ 5 . 

Alternatively, we started the field by choosing uniformly at random the phase of the order parameter between 
and 2ir, and its modulus between and p = 10 (in one set of simulations) or between and p = 10~ 6 (in another 
set). As expected, the results obtained are mostly independent of how the initial conditions are set pffl. 
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C. Numerical experiments 



The presence of topological defects (or candidate ones) in the field A(x), can be determined from the fact that for 
any closed curve C we have 

j> dx.V0(x) = 2im c (21) 
c 

where A(x) = p(x) exp[i#(x)] and nc is the total topological charge of the defects inside the curve. By candidate 
topological defects we mean those who have a net circulation of the phase, but not the equilibrium profile. That is, 
we have the phase defect, but the modulus is still evolving. 

The presence of a vortex can be observed trough expression (|2l| ) . Numerically, we will sum the shortest difference 
of the phase of the order parameter field along the lines between nodes of the grid surrounding each plaquette. Let 
us call 

s(a, (3) = (3 — a 
if(s > 7r) s = s — 2ir 
if(s < -7r) s = s + 2-k (22) 

So, four neighbor sites k, I; oriented counter clockwise will yield 

v = ^-(s(6i,6 i ) + s(e k ,6 j ) + s(6 l ,6 k ) + s(6 i ,6 1 )) =±1,0 (23) 

This device can only measure vorticity \v\ < 1. So, we can not detect a pair vortex-antivortex laying in a single 
plaquette, nor vortexes with greater vorticity. But this is enough, given the mutual annihilation of very close vortexes 
and the almost absolute absence of higher vorticity, which can be seen in the representation of the phase of the order 
parameter field. By the way, periodic boundary conditions provide a test of the accuracy of the observation, because 
the net vorticity must vanish. Higher sensitivity devices can be implemented considering higher plaquettes, which 
means more surrounding sites. 

The reciprocal representation of the field is obtained trough the usual fast Fourier transform (FFT) p§). In two 
dimensions this gives a discrete representation of A (k) at sites k = (n,m)j^ with n,m — •••) T — 1> an d a is 
the net parameter supposed equal in both directions. With our discretization, k = (n,m)j^, and the domain of the 
reciprocal representation embodies the circle k < as well as a number of higher modes. 

The various mean values of the field can be obtained easily in each time step, since the space and ensemble average 
commute. On the other hand, the power spectrum requires saving each run for further processing. 

In order to try to measure the correlation length, we consider the k 2 dependence of the ensemble dispersion g 2 of 
the amplitude of the kth mode. 

A- 



A(k) (24) 



For a thermal distribution this is a straight line, 



\ = k i±n (25) 

We can estimate — > when k 2 — > 0, through a linear fit of the ensemble media of the spectra. This is a rough 
estimate, but gives a qualitative description of the behavior of the correlation length. 

A better determination of the correlation length can be obtained from the out of equilibrium distribution of modes 
f[9| (see below). At long wavelengths k 2 < £ -2 , 

g 2 a he-^ (26) 
This factor can also be measured from a linear fit, and the system quickly reaches this regime. 
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IV. RESULTS 



A. Anatomy of a quench 

Figure [l] shows a typical evolution of a quench. We have plotted the ensemble average of the absolute value of the 
order parameter field (|A|), the HML prediction n t for the defect density, the magnitude 1/ (47r£ 2 ) (where £ is the 
correlation length measured form a fit of the long wavelength part of the correlation function, as in eq. (|2^)) and the 
observed defect density n Q , all as functions of time. This run corresponds to /i = 10~ 4 and T — 2. We have chosen 
t = as the point where the second derivative of the order parameter changes sign. 

The Figure ^ shows (|A|) , n t and n Q for all the six ensembles tested. The time scales are shifted in order to make 
all the inflection points of (|A|) to coincide. We can see that the behavior of each ensemble is essentially the same. 

Once the simulation begins, the gauge field (initially null) adjusts itself in order to follow the order parameter field, 
reacting back on it. This is what can be expected from the very different time scales involved in the equations (|l9| ) 
and (p0|), and it appears in the graph as the initial decay of the order parameter field. The observed discrepancy 
between predicted and observed defect densities is quickly smeared out by the evolution of the field. Actually neither 
of them is reliable this early in the simulation, since the prescribed thermal distribution has too much power at 
short wavelengths, which can be eliminated with a cutoff. Furthermore the algorithm to identify defects is blind to 
higher vorticity monopoles. Both errors are quickly self-corrected, though, as the density of defects decays and the 
distribution of modes becomes a gaussian function, as in eq. ( |2^ ) . This gaussianity should not be confused with that 
of the ensemble, but just refers to the shape of g(k) 2 . 

The first stage of the quench (once the transient is over) is characterized by the exponential grow of the order 
parameter field, which can be parametrized empirically as (|A|) = 2e al59t ~ 1 » 2e^ : ~ 1 (the dashed line in Figure 
H), where t is the synchronized time. While this is to be expected from the growth of the spinodal instability, it 
must be observed that we are already beyond the linear regime at this stage. This is the regime where the HML 
prediction is essentially exact. The ensemble probability density of the field is clearly gaussian, as can be tested by 

the ratio (|A| 4 \ (|A| 2 \ . In fact, for a gaussian ensemble, considering only diagonal terms of the correlation we 



have \|A| 4 y = 2 ^|A| 2 y . We have plotted this ratio vs time, in Figure |[ together with the ratio n t /n and (|A|) . 
The plot represents the average of these quantities over all six ensembles; the dot lines around the first two represent 
the dispersion between ensembles. The dashed lines around the plot of the order parameter represent the empirical 
fit to an exponential, and the tangent at the inflection point. 

When the growth of the average order parameter ceases to be exponential (around t ~ —50, see fig. Q), we enter a 
transition regime where first the gaussianity of the ensemble, and then the HML prediction cease to hold 

In the final stage, from t ~ 25 on, the approach to the equilibrium value is also exponential, as can be appreciated 
in Figure [|, where we have plotted the time derivative of (|A|), with linear-logarithmic scales. In this final stage, 
the exponential behavior of the field is oc 1 — e — 0.134^ anc j topological defects have attained almost their stable 
profiles (see below). 

We can see that, as long as the ensemble is Gaussian, the HML prediction is exact for all practical purposes. 
Around t ~ —30, both gaussianity and the exponential growth of the order parameter break down; however, the HML 
estimate is still a good approximation until later times, t ~ —10. 



B. Evolution of the structure function 



The two main regimes in the evolution of the quench, the early one dominated by the growth of the order parameter, 
and the late one dominated by the evolution of the defect network, are also clearly seen in the evolution of the 
structure function, namely the Fourier transform of the equal time order parameter correlation function (for the 
structure function in systems with global symmetry, see ref. p9[|). 

At early times and long wavelengths, the order parameter and the gauge field are essentially decoupled. Under this 
approximation, the dynamical equation eq. ([12]) becomes, 

|A=i[V 2 A + A] (27) 

Assuming that each mode has a random initial phase, the theoretical prediction for the structure function at early 
times and long wavelengths is 
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g\ ~ exp I (I j [1 - fc 2 ] J (28) 

In particular, the correlation length, determined from the scaling condition g\ ~ / grows as In Figure |^ 
we have plotted the correlation length squared as a function of time for each ensemble; the result clearly agree with 
expectations. In this regime the calculation (^) reduces to n t rs j^/tz, which, as we have seen, agrees very well with 
the observed density. 

For later times and wavelengths shorter than the average defect separation, the structure function is dominated by 
the profile of an isolated defect. The Abrikosov-Gorkov vortex centrated at the origin is given by 

A~(l-e"~^)e ie (29) 

where n is the characteristic size of the defect. The Fourier transform of this shape gives 



~ 2tt 

Afc ~F 



1- (l + rr 2 fc~ 2 ) 2 e- 4kX0 (30) 



-ik-xn 



With xo the position of the defect. Considering short wavelengths, as compared to the average defect separation, 
we obtain 



9l 



A, 



(31) 



where n a is the observed density of defects. By construction, a grid can not support singularities like that at the 
origin in (p9j), but the power law characteristic of (jjn]) is the kind of spectrum we hope to find in the final regime of 
the quench. 

Figure ^ displays the evolution of the structure function (plotted every 10 time units) as a function of ^q/c 2 , for 
random initial conditions and p = 10~ 6 ; the vertical scale is logarithmic. The early plots clearly display the gaussian 
behavior predicted by eq. (^8|). After t — —50 (bold line) the short wavelength modes begin to grow beyond the early 
times prediction. The second bold line represents the structure function at t ~ 0, that is, the beginning of the late 
times regime. The insert shows the same plot, for a wider range in wavenumber. Figure [?] shows the same for the 
thermal initial conditions, and fj, = 10~ 4 . The same behavior is obtained, being remarkable that the value of g\ for 
fc£o = 1 remains constant over the early times regime. 

In Figure ^ we have contrasted three of the structure functions shown in Figure ^ (corresponding to times t = 50, 
70 and 340) with the structure function of an isolated defect, as given by eq s. ( p0[ ) and (|3l|), given by the solid line. 
For visual effect, we have overestimated slightly the value of n in eq. (Blj) . We have set r\ = 1, as predicted by 
theory. 



C. The epoch of defect formation 

Candidate defects exist from the very beginning of the quench, as can be detected from the phase of the order 
parameter field, and predicted by the HML formula. It is the development of a well defined vortex and the pattern of 
supercurrents around it which makes the system leave the gaussian distribution. The exponential growth slows down 
as soon as the modes stop behaving almost independently (t ~ —30). The formation of the shape of the vortex itself 
is also starting at this time. 

We may follow the formation of the vortexes through the evolution of the " kinetic " free energy J30| 

K = yJ d2x IHV-A)A| 2 =1 Jd 2 x {|Vp| 2 +p 2 |V^- A| 2 } (32) 

where A = pe lS , and the last term corresponds to the supercurrents. Initially K is very low, and starts building up 
with the steeping of the field gradients around the candidate defects. When the defect attains its final shape, both Vp 
and the current die out outside the core, but there is a core contribution left, and so K reaches a final value which is 
proportional to the defect density. The subsequent evolution of K simply follows the slow decay of the defect density 
due to defect - defect annihilation. 

The Figure || shows the ensemble average of K/n (also averaged over the six ensembles considered) as a function 
of time. For comparison purposes we also show 7(|A|) (dashed line), where 7 = 4.82 is the asymptotic value of 
K/n . At time t ~ 70 the final shape and current has been reached, and the tiny fluctuations are due to transients 
corresponding to vortex annihilation. 

A remarkable implication of Figure]^ is that defect formation occurs at a relatively well defined epoch, from t ~ — 30 
to 40. 
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D. Final Remarks 



This paper attempts to answer the question of how reliable are estimates of the defect density based on the 
approximation of gaussian ensembles as applied to nonlinear phase transitions. To this effect, we have measured 
independently the gaussian prediction and the actual defect density as a function of time after an instantaneous 
quench to zero temperature in a two dimensional superconductor. The evolution of the quench goes through three 
stages, an initial one dominated by the unfolding of the spinodal instability, a final one dominated by defect - defect 
interactions, and a transitional stage when most defects are actually formed. 

We find that the gaussian estimate is essentially exact over the early regime, and continues to be accurate almost 
to the end of the epoch of defect formation. This moment in time is marked by the jump in the "kinetic "energy K 
(see eq. (|§2|)); it is also the time when the order parameter reaches about half of its equilibrium value. While the 
gaussian estimate itself is not reliable beyond this point, it is still a valid tool to fix the initial conditions of the defect 
network, whose subsequent evolution must be investigated by other means (like those in this paper, or in ref. ||). 

These results confirm theoretical expectations, but it is nevertheless satisfactory to have solid numerical proof 
of formerly theoretical conjectures. The very detailed view of the process of defect formation which is afforded by 
our simulations should also be valuable in investigating more subtle processes, such as preheating during the non 
equilibrium phase transition |31| or instabilities due to strong field effects |}2|. Also, by doing a more complete 
simulation, where we could also control the quench rate, it ought to be possible to investigate Zurek's conjecture 
about the scaling of the defect density with the quench rate |l^,[ll],[l4| . Finally, it is of interest to perform simulations 
in regions of parameter space approaching actual experimental contexts. We continue our research in these manyfold 
directions. 
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FIG. 1. Ensemble average of the absolute value of the order parameter field (|A|), the HML prediction nj for the defect 
density, the magnitude 1/ (4-7r£ 2 ]_( where £ is the correlation length measured form a fit of the long wavelength part of the 
correlation function, as in eq. (|2q) ) and the observed defect density n , all as functions of time. This run corresponds to 
H = 10 -4 and T = 2.. Initially £ differs from both the predicted and observed densities, but the rise of the gauge field smears 
out this difference. While the graph of (|A|) changes concavity at t = 0, this is not easily appreciated due to the distortion 
caused by the linear-logarithmic scales. The same curve is plotted in linear-linear scales in fig. |s| 
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FIG. 3. (|A| 4 ^) (|A| 2 ^ , as a function of time, together with the ratio n t /n and (|A|) . The plot represents the average of 
these quantities over all six ensembles; the dot lines around the first two represent the dispersion between ensembles. The dashed 
lines around the plot of the order parameter represent the empirical fit to an exponential, and the tangent at the inflection 
point. The ratios show that the predicted and observed defect densities agree very well even when the field distribution ceases 
to be gaussian. 
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FIG. 6. Power spectra measured for the ensemble p — 10 6 at equally spaced different times, as a function of £ 2 k 2 , for 
random initial conditions and p = 1CP 6 ; the vertical scale is logarithmic. The first bold curve corresponds to the time at which 
the spectrum leaves the behavior g\ ~ he~^ ^ k . The second corresponds to the inflection time. The insert shows the 
growth of short wavelength modes needed for the final shape of the defects. 
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FIG. 9. Ensemble average of K/n (also averaged over the six ensembles considered) as a function of time. For comparison 
purposes we also show 7 (|A|) (dashed line), where 7 = 4.82 is the asymptotic value of K/n , showing the final between the 
kinetic term and the defect density. 
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